Experimental signature of programmable quantum annealing 
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Quantum annealing is a general strategy for solving difficult optimization problems with the aid of 
quantum adiabatic evolution [3, 01 ■ Both analytical and numerical evidence suggests that under ide- 
alized, closed system conditions, quantum annealing can outperform classical thermalization-based 
algorithms such as simulated annealing 0, 01 • Do engineered quantum annealing devices effectively 
perform classical thermalization when coupled to a decohering thermal environment? To address 
this we establish, using superconducting flux qubits with programmable spin-spin couplings, an 
experimental signature which is consistent with quantum annealing, and at the same time inconsis- 
tent with classical thermalization, in spite of a decoherence timescale which is orders of magnitude 
shorter than the adiabatic evolution time. This suggests that programmable quantum devices, 
scalable with current superconducting technology, implement quantum annealing with a surprising 
robustness against noise and imperfections. 
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Many optimization problems can be naturally ex- 
pressed as the NP-hard problem of finding the ground 
state, or minimum energy configuration, of an Ising spin 
glass model 01, 
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where the parameters hj and Jjk are, respectively, local 
fields and couplings. The operators cr| are Pauli matrices 
which assign values {±1} to spin values {t,i}- Two al- 
gorithmic approaches designed to address this family of 
problems are directly inspired by different physical pro- 
cesses: classical simulated annealing (SA), and quantum 
annealing (QA). 

SA probabilistically explores the spin configuration 
space by taking into account the relative configuration 
energies and a time-dependent (fictitious) temperature. 
The initial temperature is high relative to the system en- 
ergy scale, to induce thermal fluctuations which prevent 
the system from getting trapped in local minima. As the 
temperature is lowered, the simulation is driven towards 
optimal solutions, represented by the global minima of 
the energy function. 

In QA P, Q the dynamics are driven by quantum, 
rather than thermal fluctuations. A system implement- 
ing QA (sl-fiot is described, at the beginning of a compu- 
tation, by a transverse magnetic field 
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The system is initialized, at low temperature, in the 
ground state of -fftrans, an equal superposition of all 2^ 
computational basis states, the quantum analog of the 



initial high-temperature classical state. The final Hamil- 
tonian of the computation is the function to be mini- 
mized, -ffising- During the computation, the Hamiltonian 
is evolved smoothly from iJtrans to iJjsing, 

H{t) = A(i)i/trans + B{t)Hisine, t € [0, T] , (3) 

where the "annealing schedule" satisfies A{0), B{T) > 
and A{T) = B{Q) = 0. If the change is sufficiently slow, 
the adiabatic theorem of quantum mechanics predicts 
that the system will remain in its ground state, and an 
optimal solution is obtained (ill [l2j . Similar transforma- 
tions with more general Hamiltonians are equivalent in 
computational power (up to polynomial overhead) to the 
standard circuit model of quantum computation [l3l . [l3 | , 
and offer at least a quadratic speed-up over any classical 
SA algorithm fisj . 

Realistically, one should include the effects of cou- 
pling to a thermal environment, i.e., consider open sys- 
tem quantum adiabatic evolution |16l42ll |. An imple- 
mentation of open system QA has recently been re- 
ported in a pro gram mable architecture of superconduct- 
ing flux qubits |22H25j . and applied to relativel y si rnple 
protein folding and number theory problems [261 . |27| . 
Althoug h q uantum tunneling has already been demon- 
strated (25[, the decoherence time in this architecture 
can be three orders of magnitude faster than the compu- 
tational timescale, due in part to the constraints imposed 
by the scalable design. In the circuit model of quantum 
computation this relatively short decoherence time would 
imply, without quantum error correction [28l[29j . that the 
system dynamics can be described by classical laws [s^] ■ 
In the context of open system QA, this might lead one to 
believe that the experimental results should be explained 
by classical thermalization, or that in essence QA has ef- 
fectively degraded into SA. 

Here we address precisely this question: are the dy- 
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FIG. 1: Connectivity graph of the degenerate Ising Hamilto- 
nian used in our experiments. The four spins in the central 
square are the "core spins" [the first four in Eq. Q], the four 
peripheral spins are "ancillae spins" [the last four in Eq. Q] . 
As depicted, the local fields hj of the core spins have value 
+1, the local fields of the ancillae spins have value —1, and 
all couplings Jjk are ferromagnetic with value 1. 



namics in open system QA dominated by classical ther- 
malization with respect to the final Hamiltonian, as in 
SA, or by the energy spectrum of the time-dependent 
quantum Hamiltonian? We answer this by studying 
an cight-qubit Hamiltonian representing a simple opti- 
mization problem, and show that classical thermalization 
and QA make opposite predictions about the final mea- 
surement statistics. Our Ising Hamiltonian, depicted in 
Fig. [2 has a 17-fold degenerate ground state 

{ itttn;;;) , • • • , itnnn;) , • • • , inntnt) } (4a) 
umm) , (4b) 

Sixteen of these states form a cluster of solutions con- 
nected by single spin-flips of the ancillae spins [Eq. (|4a|) ]. 
while the seventeenth ground state is isolated from this 
cluster in the sense that it can be reached only after at 
least four spin- flips of the core spins [Eq. (j4bp ]. As we 
show below, classical thermalization predicts that the iso- 
lated solution will be found with higher probability than 
any of the cluster solutions, i.e., it is enhanced. Fur- 
thermore, after an initial transient, faster thermalization 
corresponds to a higher probability of flnding the isolated 
solution. Open system QA makes the exact opposite pre- 
diction: after an initial transient, the isolated solution is 
suppressed relative to the cluster, and faster quantum 
dynamics yields higher suppression (lower probability). 
Our experimental results are consistent with the open 
system QA prediction of the suppression effect, and in- 
consistent with classical thermalization. We next discuss 
these opposite effects, starting from the classical case. 

Let Pi denote the probability of state i in the clus- 
ter (|4ap . and p^ the probability of the isolated state (|4bp . 
The thermalization dynamics are dominated by single 
spin- flips in our experiment (see Appendix). The prob- 



abilities Pi are all close because states in the cluster are 
connected by single spin flips^ so we consider the average 
cluster probability pc = X]i=i Pi/16. Enhancement of 
the isolated state means that ps > pc- Our SA numer- 
ics show that this is indeed the case for different update 
rules and cooling schedules, throughout the thermaliza- 
tion evolution (see Appendix). To explain this, note that 
the general features of a thermalization process are de- 
termined by the spectrum of i?ising and by the combi- 
natorics of state interconversion. Each of the 17 degen- 
erate ground states can be reached from any other state 
without ever raising the energy via a sequence of single 
spin-flips, so that SA never gets trapped in local minima 
(Appendix) . 

The SA master equation and the classical thermaliza- 
tion prediction ps > pc can be derived from first prin- 
ciples from an adiabatic quantum master equation [21 1. 
Let Hs{t) and Hsb = J2a ® denote the system 
and system-bath Hamiltonians. The Lindblad equation 
is 

p=-i[Hs,p] (5) 

a/3 a^b ^ 

+ X! X! ^"-S iaa,/3piL,a " \ { L\a,aLbb,0 , p] 

aj3 ab 

where Lab,a = \a){a\Aa\b){b\, ujab = Eb - Ea, {\a)} is 
the instantaneous eigenbasis of Hs for spin vector a, and 
= JZ^dTe'''^{Bl{T)Bf3{0)). We are interested 
in the thermalization process in which the density oper- 
ator is diagonal in the computational basis of spin vec- 
tors. The system-bath coupling Hamiltonian then has 



the form Hsb = T,r€{+-,z}T,f=i9f''^j <^ -^f^' '^'^'^^'^ 
a"^ ~ {a^ ± iay)/2. We denote by (a~) the spin vec- 
tor resulting from flipping the j th spin up (down) . From 
here we arrive at the classical master equation for the 

populations Pa = Paa- 

N 

j=l r=± 

(6) 

and the detailed balance condition f{Ea — E^±) = 

exp[-/3(S^± - Ea)]f{E^± - Ea). Eq. ^ is the' mas- 

ter equation that we used in our SA numerics. It can 
also be used to derive the classical thermalization pre- 
diction Ps > PC- To this end, it can be seen directly 
that the isolated state is connected via single spin-flips 
to 8 excited states with energy —4, giving the rate equa- 
tion Ps = 8/(— 4)pe — 8/(4)ps. In contrast, all states in 
the cluster are connected via single spin-flips to at most 
4 singly-excited states; the other four spin-flips connect 
between other states in the cluster and hence conserve 
the energy. Thus pc ~ 2 {f{-A)pe - f{'i)pc)- Compar- 
ing, we conclude that population feeds into the isolated 
state faster than into the cluster and, given that initially 
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FIG. 2: The time-dependent gap between the ground state 
and the lowest six excited states in the relevant region of 
the experimental QA evolution. After time t = 0.5T the 
highest energy level shown corresponds to the isolated state. 
The inset shows the transverse field magnitude A{t) and Ising 
Ifamiltonian magnitude B{t) used in our experiments, during 
the same time interval. 



Ps ^ PC, we always observe Ps >Pc- Simultaneous dou- 
ble spin-flips do not change this conclusion, and higher 
order simultaneous spin-flips are physically less likely. A 
complete derivation is given in the Appendix. 

We next analyze the corresponding predictions of QA. 
A crucial difference respect to SA is that now the 
relevant energy spectrum is given by a combination of the 
final Ising Hamiltonian and the transverse field. Conse- 
quently, as shown in Fig. [21 the degeneracy of the ground 
space is lifted for times t < T. The isolated state has sup- 
port only on the highest eigenstate plotted during the sec- 
ond half of the evolution. Given that the system starts in 
the ground state, the isolated state is suppressed by the 
energy gap, until this gap vanishes at the end of the evo- 
lution. The isolated state remains suppressed nonethe- 
less, since transitions to other low energy states require 
at least four spins-flips. The transverse field term, which 
drives simultaneous spin-flip transitions, is small at large 
t. If the four spins-flips arc not simultaneous, these tran- 
sitions involve excited states with much higher energy, 
and are suppressed. This predicted QA suppression of 
the isolated state is confirmed by our closed and open 
system quantum dynamical simulations. 

Our experiments were performed using the D-Wave 
One Rainier chip at the USC Information Sciences In- 
stitute, comprising 16 unit cells of 8 superconducting 
flux qubits each, with a total of 108 functional qubits. 
The couplings are programmable superconducting induc- 
tances. The qubits and unit cell, readout, and control 
have been described in detail elsewhere [22H25j . The ini- 
tial energy scale for the transverse field is lOGHz (the 
A function in Fig ^ , the final energy scale for the Ising 
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FIG. 3: Left: schematic depicting the maximal connectivity 
graph (-^"4,4) of the qubits inside a unit cell. Right: an em- 
bedding of i/ising from Fig. [T] (right). 



Hamiltonian (the B function) is 5.3GHz, about fifteen 
times the experimental temperature of 17mK w O.SSGHz. 
To gather our data, we ran each of the 144 embeddings 
4000 times, in batches of 1000 readouts, resetting all the 
local fields and couplers after each batch. 

A diagram of the experimentally achievable coupling 
configurations is shown in Fig. [3] (left). The experimen- 
tal results arc shown in Fig. |4l The key finding that 
is immediately apparent is that the isolated state is ro- 
bustly suppressed, in agreement with the QA but not the 
SA prediction. 

Is it possible that suppression has an explanation other 
than QA? The main physical argument along these lines 
is that a systematic or random bias due to experimental 
imperfections breaks the 17- fold ground state degeneracy 
and energetically disfavors the isolated state, thus lower- 
ing Ps if the system thermalizes. We proceed to examine 
this and the robustness of the suppression effect. 

First, note that spin numbers j = 1,...,8 must be 
assigned to the flux qubits before each experimental 
run. One of the 144 possible such "embeddings" al- 
lowed by the symmetries of the Hamiltonian and the 
hardware connectivity-graph is shown in Fig. [3] (right). 
Second, note that spin-inversion transformations H{t) — > 
ajH{t)aj commute with -fftians, and simply relabel the 
spectrum of both i?ising and H{t): if a certain spin con- 
figuration has energy then the corresponding spin 
configuration with the jth spin flipped has the same en- 
ergy E under crj iJisingcJ. Spin-inversions also commute 
with the spin-flip operations of classical thermalization. 
Therefore all of our arguments for the suppression of the 
isolated state in QA and for its enhancement in classi- 
cal thermalization are unchanged. Using spin inversions 
we can check that the suppression effect is not due to a 
perturbation of the Hamiltonian such as a magnetic field 
bias. Indeed, by performing a spin-inversion on all N 
spins we obtain a new Ising Hamiltonian where the iso- 
lated state is that with all spins-up. If a field bias sup- 
pressed the all spins-down state, then it would enhance 
the all spins-up state. Figure [5] rules this out. We also 
tested cases with only antiferromagnetic couplings, and 
with random spin-inversions. The results for one such 
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FIG. 4: Statistical box plot [sJl of the experimental results for ps (left columns) and pc (right columns). The total annealing 
time was T — 5/iS in each run. In each column the bar is the median, the box corresponds to the lower and upper quartiles, 
respectively, the segment contains most of the samples, and the +'s are outliers. Cells 1-4 are physically distinct 8-qubit unit 
cells on the chip. No statistically significant variation is seen as a function of the unit cell number. "Inv" is the case where each 
local field hj is fiipped to —hj. In this case the isolated state corresponds to the state [tttttttt)) the opposite of Eq. (I4b|l . 
While this has a small effect for cells 2 and 3, in all cases the isolated state is significantly suppressed, as predicted by QA. 
This establishes that suppression of the isolated state is not due to a global magnetic field bias. 



random inversion example are shown in Fig. [51 In all 
cases we found agreement with the QA prediction, but 
not with classical thermalization. 

Robust suppression holds even at the level of individual 
embeddings and spin inversions. We found that Ps < 
3%, while PC ^ 6% for each of the thousands of such 
cases we tested (the highest median for the experiment 
in Fig |4] is 0.004). Thus suppression survives breaking 
of the ground state degeneracy, which certainly occurs 
due to the limited precision of '--^ 5% in our control of 
{hj^Jjk}. The suppression effect is robust because it 
does not depend on the exact vakies of these parameters, 
but on the relatively large Hamming distance between 
the isolated state and the cluster. 

Finally, we consider the effect of increasing the anneal- 
ing time. Open quantum and classical systems converge 
towards thermal equilibrium. Therefore if the cause of 
suppression is the QA spectrum, longer annealing times 
will result in ps increasing^ approaching its Gibbs dis- 
tribution value. This would not be the case if ps were 
governed by the spectrum of -ffising- In Fig- [H we com- 
pare a numerical simulation of open system QA, using an 
adiabatic Markovian master equation [2l| , with classical 
thermalization. The quantum prediction of increasing ps 
is confirmed experimentally, as shown in Fig. [S] 

We thus arrive at our main conclusion: signatures of 
QA, as opposed to classical thermalization, persist for 
timescales that are much longer than the single-qubit 
decoherence time (from 5//s to 20ms vs tens of ns) in 
programmable devices available with present-day super- 
conducting technology. Our experimental results are also 
consistent with numerical methods that compute quan- 



tum statistics, such as Path Integral Monte Carlo ^. Our 
study focuses on demonstrating a non-classical signa- 
ture in experimental QA. Different methods are required 
to address the question of experimental computational 
speedups of open system QA relative to optimal classical 
algorithms. 
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Appendix A: First-principles derivation of the 
master equation 

Here we derive the master equation used by SA from 
first principles within the open quantum systems formal- 
ism. This motivates classical SA as a model for a system 
dominated by classical thermalization of the final Ising 
Hamiltonian. 

Let Hs{t) be the time-dependent system Hamiltonian 
and HsB = J2a ® ^^'^ system-bath Hamilto- 

nian. We have previously established that the Lindblad 
equation within the rotating wave approximation has the 



^ M. Troyer and T. Roennow, private communication. 
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FIG. 5: Statistical box plot of the probability of the isolated state for a fixed set of qubits, with different annealing times. 
In this plot each of the 144 possible embeddings is averaged with the same embedding after a complete spin inversion. This 
compensates for global magnetic field biases (which can be seen in Fig. |4l cells 2 and 3). The Ising Hamiltonian for this data 
was obtained by applying a random spin inversion to Huing from Fig. [T] The probability of the isolated state increases with the 
annealing time T, in contrast to the classical thermalization prediction. While in classical thermalization an initial distribution 
concentrated around the cluster can also result in suppression of the isolated state, this suppression is highly unlikely to persist 
after a random inversion. 



form 

p^-i[Hs,p] (Al) 
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al3 a^b 



ab 



Laa,l3pL 



6b, ( 



{Lia,aLbb,l3,p} 



where 



Lab,a = \a){a\Aa\b){b\ 
Lib^a = \b)mi\a){a\ 

^ab = Eb — Ea , 



(A2a) 
(A2b) 
(A2c) 



{|a)} is the instantaneous eigenbasis of Hs (we have 
suppressed its explicit time-dependence) for spin vector 
a = {ai, . . . , ajv}, where G {t, i} , and 



dre-Mi?t(r)B^(0)) (A3) 



is the Fourier transform of the bath correlation function. 
The star adornment on the first subscript {a*) is a re- 
minder that the first operator in the bath correlation 
function is Hermitian-transposed. We have ignored the 
Lamb shift in Eq. since for a time-dependent Lind- 
blad evolution it amounts to a small perturbation of the 
system Hamiltonian. We used this form of the master 
equation for our quantum open system numerical simu- 
lations, as detailed elsewhere I21j. 



We show in Sec. [B]that for a bath in thermal equilib- 
rium at inverse temperature (3 



7a*/3(-w) = e P"7/3Q*(a;) , 



(A4) 



where 



dre--(B^(r)i?t(0)) . (A5) 



We assume that the system-bath coupling Hamiltonian 
has the form 



N 



j=i re{±.z} 



(A6) 



where cr^ = (CT^±iCT^)/2, we identify |t) with |0) and ||) 
with and where we neglect higher-order interactions 

of the form crj ® ct^ (gi -Sj^*' or above. Since Hsb is 

Hermitian we also have B^^^^ = Bf^\ s'f'^'^ = B'f\ 

dj"^^* ~ sj^'' 9j^^* — 9j^^^ ^'^'^ where the asterisk denotes 
complex conjugation. In the computational basis of spin 
vectors {a}, we introduce the notation 



(A7) 



which denotes either a flipping of aj, or if either 
acts on aj or aJ' acts on aj =^. Then 

(-H\b)-{-f).b-Kbf^ (A8) 
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FIG. 6: Probability of the isolated state for numerical simula- 
tions of classical thermalization (Metropolis update rule) and 
open system QA as a function of the total annealing time T. 
"Ideal" vs "perturbed" corresponds to simulations for i/ising 
without and with a perturbation which increases the energy of 
the isolated state (red). In classical thermalization ps always 
decreases with T, while it increases for QA in the ideal case. 
It remains almost constant for QA with the perturbed Hamil- 
tonian. Even if the isolated state is suppressed energetically 
due to a perturbation of //ising, fast classical thermalization 
can still enhance its probability. QA with the ideal Hamilto- 
nian gives the best qualitative fit to the experimental data. 
System-bath coupling in the QA simulation corresponds to a 
decoherence time of 150ns. 



where the 6 function is defined to evaluate to zero also 
when annihilates \b). We are interested in classical 
thermalization, in which the density operator is diagonal 
in the computational basis {a}, so we set pat ~ for 
a ^ b. Equation (|A1[) then gives pab = 0. Using indexes 
c = {f,j) and /? = (s, k) in Eq. (|A1[) . where r, s € {±, z} 
and j, fc G [1,...,A^], and taking the diagonal (a| • \a) 
matrix element, the Lindblad equation becomes 



E l{r.jY{s.k){^ab) K)abPbb [i'^j )^) ba 
b\b^a 



vanish, and using Eq. (jASp . we are left with 

N 

P, 



l(r,jy(r,j){^a-a)Vc^ , (AIO) 



where we denoted Pa = Paa-, the probability of spin con- 
figuration a. We can furthermore identify 



(Alia) 
(Allb) 



as the transition probabilities, so that Eq. (jAlOp becomes 
the rate equation 

N 

P« = E E P^^7 ^ «K7'- ~ ^ a^,)Va . (A12) 

j=\ r=± 

This can be further simplified using the KMS con- 
dition. Indeed, note that, using Bair) = '^'f^'^) 
Eqs. dMI and ((X5|) . we have 

7(±jr(±.j)(^) = 7(T.j)(Tjr(^)- (A13) 

Using this along with w^±^ = — w^^± [Eq. (jA2cP ] and 

i-i— t ^ ^ 

Eq. (|A4p . we have 

-(3w ± 

7(±j)*(±j)K±a) =e 7(Tj)*(Tj)Ka±) • (^14) 

Therefore Eq. (jAlip yields 

^ 4) = 7^^^^^ Isf ^P7(T,.r (T,.)Ka±) (Al5a) 
P(a± ^ a) = Isf 'l'7(T..)*(T,.)Ka±) • (A15b) 

This, together with .9^'''''* = gives the detailed bal- 
ance condition for thermalization dynamics 



P(a ^ a± ) - p (is^± /j {^a - E^± ) 



, (A16) 



where we introduced transition functions fj{AE), which 
we identify with the transition probabilities in Eq. (jA15|l . 

We can now rewrite Eq. (|A12p as the classical master 
equation that we used in our SA numerical simulations 



N 

Pa = E E - - /^ (^- - 



j=l r=± 



(A17) 



Note that the sum in Eq. (|Aip involving the resonant con- 
tribution 7q*^(0) vanishes, since the terms Laa,[jpL\j^ ^ 

and ^ {^QQ a-^f>fc,/3i P} cancel after taking the diagonal 
matrix clement. Moreover, since Eq. (|A9p involves only 
off-diagonal terms (6 7^ a), contributions due to cr^ all 



Appendix B: Correlation functions and the KMS 
condition 

Here we derive the detailed balance condition Eq. (jA4p 
from first principles. Our calculation closely follows 
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Ref. |2l|, but differs in that it applies also to non- 
Hermitian bath operators. 

The correlation function of a thermal bath is as- 
sumed to satisfy the KMS (Kubo-Martin-Schwinger) con- 
dition ii 



(i?t(r)i?^(0)) = (i?^(0)i3t(r + zP)) 



(Bl) 



This expression has the advantage that it also applies 
to operators which are not trace class. For trace class 
operators the KMS condition can be derived assuming 
that the bath is in a thermal state, pB = e~P^^, where 
Hb is the bath Hamiltonian. In this case: 

{Bi{T)Bpm = Tv[pBui{T,o)BiuB{T,mp] 

= lTr[B^,e-(P-*^)^«Bte-*^^«] 

= i-Tr[B,3e'^'"+''^'^^B^e"''^^+*P'-^^e"''-^^] 



= Tr[pBBpUl{T + z|3, 0)St C/b(t + i|3, 0)] 
-(i?M0)i?t(r + *|3)) , 



(B2) 



where J7b is the bath unitary evolution operator. Note 
that 



{Bi{T)Bpm = {Bp{-T-zmim ■ 



(B3) 



If in addition the correlation function is analytic in the 
strip between r = — i|3 and t = 0, then it follows that 
the Fourier transform of the bath correlation function 
satisfies the detailed balance condition Eq. (|A4p as we 
show next. 




FIG. 7: Contour used in our proof of the KMS condition. 
We compute the Fourier transform: 

/•OO 

7o*/3(^) - / dre--(Si(r)B^(0)) 

) 

dre-^(S0(-T- iP)i3t(O)) . (B4) 



To perform this integral we replace it with a contour inte- 
gral in the complex plane, §^ dTe*'^'^(i3^(— r — i(3)i?J(0)), 
with the contour C as shown in Fig.[71 This contour inte- 
gral vanishes by the Cauchy-Goursat theorem [ssj since 



the closed contour encloses no poles (the correlation func- 
tion {Bp{t)BI^{0)) is analytic in the open strip (0 —i(3) 
and is continuous at the boundary of the strip [3J]), so 
that 



(...) = 



(B5) 



(•••)+ / (•■•)+ / (•■■)+ / (• 

t Ji 



where (...) is the integrand of Eq. (|B4p . and the integral 
is the same as in Eq. (|B4p . After making the variable 
transformation r = —a; — «|3, where x is real, we have 



(...) = -e 



e-^-^B0{x)Bim . (B6) 



Assuming that {Ba{ioo — i^)Bp{0)) = (i.e., the corre- 
lation function vanishes at infinite time), we further have 
(...) = J^{. . .) ~ 0, and hence we find the result: 



dTe--(S^(-r-^P)St(0)) (B7) 

OO 

eP- / e— (B^(r)i3t (0)) = e^-jp^.{-oj) , 



which, together with Eq. (|B4|, proves Eq. (|A4|) . 



Appendix C: Spectrum and ground states of the 
Ising Hamiltonian 

The spectrum of the 8-qubit Ising Hamiltonian we con- 
sider in the main text. 




can be analyzed by first considering the spectrum of the 
Hamiltonian coupling a single ancilla spin to a core spin, 
i.e., 



+1 



-1 



+1 
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The spectrum, with the core (ancilla) spin written first 
(second) is 



Itt) 

in) 
i;t) 

lU) 



-1 
-1 

3 

-1 



(CI) 



The minimum energy is —1 whether the core spin is up 
or down. It is important to note that if the core spin is 
up, the minimum energy is —1 whether the ancilla is up 
or down; this will give rise to a 16-fold degeneracy when 
we account for all spins below. 

We analyze the core spins' energies by first taking into 
account only their couplings. That is, we analyze the 
ferromagnetic Hamiltonian 



V 



Denoting by s the number of satisfied couplings (both 
spins linked by the coupling have the same sign) , the en- 
ergy is 4 — 2s, where s £ {0,2,4}. The ground states 
of this Hamiltonian are the configurations | tttt) and 
Since Eq. (|C1|) shows that the minimum energy 
of a corc-ancilla pair is —1, when adding the low energy 
configurations of the couplings to the ancillae the mini- 
mum energy is —8. It also follows from Eq. (jCip that the 
ground state configurations of the full 8-qubit Hamilto- 
nian are 



mmm , 



(C2a) 
(C2b) 



where the first (last) four spins are the core (ancillae) 
spins, and |^) means that the spin can be cither up or 
down. The all-spins down case (|C2a[) results from the HI) 
configuration in Eq. (|Cip . while the 16-fold degenerate 
case (IC2bp results from the degeneracy of |tt) and||t)- 

An important feature of the energy landscape of the 
8-qubit Hamiltonian is that it does not have any local 
minima. This can be easily proved by showing that a 
global minimum can always be reached from any state 
by performing a sequence of single spin flips and never 
raising the energy. To see this, consider an arbitrary state 
of the system. We can first flip all the ancillae spins to 
ID which, according to Eq. (jCl|) . can be done without 
raising the energy (independently of the state of the cor- 
responding core spin) . Then we can flip the core spins in 
order to satisfy all the couplings between core spins, ei- 
ther making them all |^) or all| f) , whichever requires the 
fewest spin flips. Again, according to Eq. (|C1|) . this op- 
eration will not raise the energy of the core-ancilla pair. 
Hence, the flnal state is either the isolated ground state 
mmm), or the state Ittttiiii) that belongs to the 
degenerate cluster of ground state configurations. 



Appendix D: Simulated annealing and classical 
thermalization 



Here we report the results of our numerical simula- 
tions of classical thermalization (or SA), using the master 
equation (|A17|1 . In the plots below we used the Metropo- 



lis update rule for the transition probability P{a 



Explicitly, if AE = E ± — Ea is the energy difference for 
the update, the transition probability is 



# of spins 



min(l,exp(-^Ai;)) 



(Dl) 



We have also tested other update rules, such as 
Glauber's [sll, and the results are essentially unchanged. 
The result are also essentially unchanged when using a 
transition probability pmin (1, exp(— /3A£')) for positive 
p (such as when simulating a continuous time master 
equation). We do find a dependence on the choice of 
annealing schedule, i.e., the functional dependence of the 
temperature on the number of steps. Three different an- 
nealing schedules wc used are shown in Fig. [8l and the 
corresponding SA results are shown in Fig. [SI The prob- 
ability ps of the isolated state is always above the average 
probability pc for a state in the cluster. 

It might be argued that thermalization at constant 
temperature corresponds most closely to the experimen- 
tal situation, given that the experimental system remains 
at an almost constant 17mK. This is modeled by the ex- 
ponential annealing schedule, which rapidly converges to 
a nearly constant temperature, as can be seen in Fig. [8] 
On the other hand, the energy scale of the Ising model 
changes during the QA evolution (see the Fig. 2 insert 
in the main text), and the cooling schedule is deter- 
mined not by the temperature alone but rather by the 
ratio between the energy scale and the temperature. We 
also show Ps and Pc for an exponential schedule with 
Metropolis updates and different numbers of steps in 
Fig. [101 



Appendix E: Classical master equation explanation 
for the enhancement of the isolated state 

We now explain why, as seen in the numerical simu- 
lations shown in Figs. [9] and [TOl the probability of the 
isolated state never exceeds that of the average of the 16 
cluster ground states, i.e., why 



1 ^ 

^16 



(El) 



We are interested in sufficiently slow thermalization pro- 
cesses (relative to spin flip rates), so that states connected 
by single spin-flips have similar populations. The cluster 
of 16 degenerate ground states and the isolated ground 
state are connected via a plateau of excited states with 
energy —4. 



9 



6 -> 
0) I \ 
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- Ti/exp(step) 

- Ti/log(step) 
Ti/step 



200 300 
Annealing step 



FIG. 8: Temperature as a function of annealing step n for 
three different schedules: exponential T{n) = Tir"^p, linear 
T(n) = Ti/{nrun + 1), and logarithmic T{n) = Ti/(log(n + 
l)nog + 1), where Ti is the initial temperature, T/ is the final 
temperature, and Texp = (T/ZT,)'''"'"' , ni„ = [Tf /T,-l) /ritot 
and riog = {Tf /Ti - f)/log(ntot + 1), where ntot is the total 
number of annealing steps. 




200 300 
Annealing step 



FIG. 9: Probabilities from SA for the three different schedules 
shown in Fig. [S] The probability of the isolated state ps is 
always higher than the average cluster state probability pc- 



Let us first derive a rate equation for the isolated state. 
A single spin-flip of a eore spin in the isolated state raises 
its energy by 4, sinee it violates two couplings between 
the core spins and corresponds to a transition from |||) to 
Iti) (where the second, ancilla, spin is unchanged), which 
doesn't change the energy according to Eq. (|C1|) . Like- 
wise, a single spin-flip of an ancilla spin in the isolated 
state violates no couplings and corresponds to a tran- 
sition from lit) to lit) (with the core spin unchanged), 
which raises the energy by 4 according to Eq. (jCl[) . There 
are 8 ways this can happen (4 core and 4 ancilla spins 
can be flipped) . Since this accounts for all the single spin 
transitions, Eq. (|A17|) yields the rate equation 



Ps = 8/(-4)pe - 8/(4)p, 



(E2) 



0.14 



0.12 



0.1 - 



■ 0.08 



0.06 



0.04 - 



0.02 - 



0.014 



0.012 




5 4 3 

Temperature 



FIG. 10: Probabilities from SA for varying total numbers of 
steps. We used the Metropolis update rule with an expo- 
nential schedule. The lines correspond to 100 (dotted), 1000 
(dashed) and 10000 (solid) steps. The upper curves (blue) 
correspond to ps, while the lower curves (red) correspond to 
PC- The inset is a magnification of the boxed part. The sepa- 
ration between the probabilities of the isolated state and the 
cluster increases as the temperature decreases. 



rate is the same for all sites [corresponding to assuming 
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(r) 



in Eq. (lA6|l ]. 



where Pe is the population of the excited states with 
energy —4. Here we are assuming that the spin flip 



We next derive the rate equation for the cluster, once 
again accounting only for single spin flips. For states in 
the cluster the core spins are all up, and ancilla-spin flips 
are energy-preserving transitions between states in the 
cluster. For core-spin flips we need to analyze two dif- 
ferent situations. The first is a configuration in a ground 
state where the core-ancilla pair starts as |tt) and the 
core spin flips, so the state becomes | It)- This vio- 
lates two couplings, with energy cost 4, and according to 
Eq. (jClj) the energy difference between these two states 
is 4, so the overall result is an excited state with energy 
0. The second is a configuration in a ground state where 
the core-ancilla pair starts as | ft) and again the core 
spin fiips, so the state becomes This again violates 
two couplings, with energy cost 4, but costs no energy 
according to Eq. (jCl[) . so the overall result is an excited 
state with energy —4. 

Thus, a state with / ancillac with spin J, and 4 — I 
ancillae with spin f connects (via single spin-flips) to I 
excited states with energy —4 and 4— Z excited states with 
energy 0. To write a rate equation for pc = 
we shall assume that all excited states with energy (—4) 
have probability p{0) (pe)-, and all states in the ground 
state cluster have equal probability pc- Summing over 
the number I of ancilla with spin J, for each cluster state. 
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the rate equation is 
16 4 . . 

Ep^ = E J(^/M)p,-;/(4)pc 

i=l 1=0 ^ ^ 

+ (4 - 0/(-8) p{0) - (4 - 0/(8) Pc) (E3a) 
= 32(/(-4)p,-/(4)pc 

+ /(-8)p(0)-/(8)pc) , (E3b) 

so that 

Pc = 2(/(-4)pe-/(4)pc+/(-8)p(0)-/(8)pc) ■ (E4) 

For most temperatures of interest, relative to the en- 
ergy scale of the Ising Hamiltonian, the dominant transi- 
tions are those between the cluster and states with energy 
-4. Transitions to energy states are suppressed by the 
high energy cost, and transitions from energy states to 
the cluster are suppressed by the low occupancy of the 
energy states. Then 



pc«2/(-4)pe-2/(4)pc 



(E5) 



QA evolution as that of a closed system evolving back- 
ward in time). According to standard degenerate per- 
turbation theory, the perturbation Pg of the ground sub- 
space is given by the spectrum of the projection of the 
perturbation on the ground subspace. Denoting by 



(i;) air' + E mmm) (nmmi (fi) 



the projector on the 17-dimensional ground subspace, we 
therefore wish to understand the spectrum of the opera- 
tor 



Pn 



E-l 



(F2) 



The isolated state is unconnected via single spin flips 
to any other state in the ground subspace, so we can write 
this operator as a direct sum of acting on the isolated 
state and the projection on the space Ilg — Ilo~{\l){\.\)^^ 
of the cluster 



-0 ® n' - 



n' 



• 



(F3) 



To show that pg > pc, assume that this is indeed the 
case initially. Then, in order for pc to become larger 
than ps, they must first become equal at some inverse 
annealing temperature |3': Ps((3') = Pc(P') = Pg, and 
it suffices to check that this implies that ps grows faster 
than PC- Subtracting Eq. ((E5)) from Eq. (|E2|) yields 



Ps 



Pc = 6(/(-4)pe-/(4)pg) 



= 6/(-4)p, - 



P{g ^ e) 



Pa g) 



(E6) 



where in the second line we used Eq. ([XTH)) . Now, be- 
cause the dynamical SA process we are considering pro- 
ceeds via cooling, the ratio between the non-equilibrium 
excited state and the ground state probabilities will not 
be lower than the corresponding thermal equilibrium 
transition ratio, i.e., — > pS^^^I = e^'^^ . Therefore, as 
we set out to show, 



Ps-Pc>0 , 
implying that at all times ps > pc- 



(E7) 



While acting on any of the four ancillae connects two 
cluster ground states, cr^ acting on any core spin of a clus- 
ter state is projected away. Therefore the perturbation 
is given by the operator 



Pr, 



E-J 

i=5 



(F4) 



where the sum is over the four ancillae spins. 

Denoting the eigenbasis of cr^' by |±) = (|t) ± \i))/V^, 
with respective eigenvalues ±1, the transverse field splits 
the ground space of i?ising lowering the energy of |tttt 
+ + and the four permutations of |— ) in the ancillae 
spins of Itttt + + +—). None of these states overlaps with 
the isolated ground state, which is therefore not a ground 
state of the perturbed Hamiltonian. Furthermore, after 
the perturbation, only the sixth excited state overlaps 
with the isolated state. The isolated state becomes a 
ground state only at the very end of the evolution (with 
time going forward), when the perturbation has vanished. 



Appendix G: The quantum Singular Coupling Limit 
does not agree with the experimental results 



Appendix F: Degenerate perturbation theory 
explanation for quantum suppression of the isolated 

state 

We can understand the splitting of the degenerate 
ground subspace of the Ising Hamiltonian iJising by treat- 
ing the transverse field -fftrans = ~ Sj=i '^j ^ pertur- 
bation of the Ising Hamiltonian -ffising (thus treating the 



Interestingly, an open system QA master equation in 
the singular coupling limit (SCL) yields results in qual- 
itative agreement with classical thermalization, and op- 
posite to our weak couplin g lim it (WCL) master equation 
(|Aip . Here, following Ref. |32, we present a derivation of 
the SCL master equation. 

We consider a Hamiltonian of the form: 

Hit) = Hs{t) + e-^Hi + e-^HB , (Gl) 
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where we take the interaction Hamiltonian Hj to have 
the form A ^ B, where the system (A) and bath (B) 
operators are both Hermitian. The formal solution in 
the interaction picture generated by Hs and Hb is given 
by: 



p{t) = p(0) - ie"^ / ds 
Jo 



(G2) 



Plugging this solution back into the equation of motion 
and taking the partial trace over the bath, we obtain: 



dt 



Psit) 



dsTif 



Hi{t), Hi{s),p{s) 



(G3) 

where we have assumed that Ti[pbB] = (B) = 0. Under 
the standard Markovian assumption that p{t) = ps{t) <E) 
Pb and under a change of coordinates s = < — r, we can 
write: 



dt 



Ps{t) 



dT [{A{t)p{t - T)A{t 



(G4) 



^A{t)A{t - T)p{t - t)) {B{t)B{0)) + h.c] 



where A{t) ~ Us{t)AUl{t) and where we have used the 
homogeneity of the bath correlation function to shift its 
time-argument. We change coordinates r = e^r' and ob- 
serve that under this coordinate change {B{t)B{0)) is 
independent of e. We assume that this bath correlation 
function decays in a time tb that is sufficiently fast, such 
that Tb <C t/e^. This allows us to approximate the in- 
tegral by sending the upper limit to infinity. We also 
assume that tb <C t'e^, which forces the correlation time 
of the bath to zero, hence its spectral density to become 
flat, and hence — using the KMS condition — amounts to 
taking the infinite temperature limit. Under these as- 
sumptions, we can now take the e — > limit, yielding the 
SCL master equation 



-i[ffs(t) + FLS,pW] 

+7(0) [Ap{t)A-\{A\p[i)) 



where 



7(w) = 



dr- 



.1 ^ — iujr' 



B{T')Bm , 



-A' I dLU-f{Lj)V ( ^ ) , 



(G5) 

(G6) 
(G7) 



where iJLS is the Lamb shift (renormalization of the sys- 
tem Hamiltonian) and where V denotes the Cauchy prin- 
cipal value. Thus, even if Hs is time-dependent, we re- 
cover the same form for the SCL master equation as in 
the time- independent case [32|. This SCL master equa- 



tion yields results in qualitative agreement with classical 
thermalization, and opposite to the WCL presented in 



Fig. 6 of the main text. Namely, the SCL master equa- 
tion predicts that the isolated state is enhanced relative 
to the cluster of states. Since the bath and system-bath 
coupling dominate the system Hamiltonian in the SCL, 
it predicts that decoherence takes place not in the instan- 
taneous energy eigenbasis of the system but in the com- 
putational basis. Furthermore, the resulting Lamb shift 
in this limit, diagonal in the computational basis, pref- 
erentially lowers the energy of the all-up and all-down 
states relative to the other computational states. To- 
gether, these two effects cause the isolated state to be 
more populated than the average population of the clus- 
ter at the end of the evolution, in contradiction to our 
experimental findings. The agreement we find between 
our experimental results and the WCL master equation 
instead supports the idea that decoherence takes place in 
the instantaneous energy eigenbasis of the system and/or 
that the Lamb shift docs not dominate the system Hamil- 
tonian. 



Appendix H: Entanglement 

Another interesting question is the possible presence of 
entanglement during the evolution. Answering this ques- 
tion experimentally requires measurements to be per- 
formed during the annealing, a capability that is absent 
from the device used in our experiments. However, when 
we compute the concurrence of the states obtained from 
our master equation (|Aip during the QA evolution, which 
is consistent with the statistics of the measured output, 
we find it to be finite, as seen in Fig. [TlJ This suggests 
that entanglement is being generated in our experiments. 



External edge (thermal) 

Center edge (thermal) 

External edge (ground state) 

Center edge (ground state) 




FIG. 11: Concurrence generated during our QA simulations, 
between a pair of ancillae qubits ( "external edge" ) and a pair 
of core qubits ("center edge"). We show the concurrence for 
the ground state and for the Gibbs state ("thermal"). Our 
master equation (|AH) for the time-dependent density matrix 
gives concurrence values between these two extremes, that de- 
pend on the system-bath coupling strength used in the simu- 
lation. 
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